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f^ . Abstract 

The question of the nature of the dark matter in the Universe remains one of the 

C^h' most outstanding unsolved problems in basic science. One of the best motivated particle 

JL , physics candidates is the lightest supersymmetric particle, assumed to be the lightest 

X I ■ neutralino - a linear combination of the supersymmetric partners of the photon, the Z bo- 

"^ ' son and neutral scalar Higgs particles. Here we describe DarkSUSY, a publicly-available 

j^ , advanced numerical package for neutralino dark matter calculations. In DarkSUSY one 

can compute the neutralino density in the Universe today using precision methods which 

include resonances, pair production thresholds and coannihilations. Masses and mix- 

^^ , ings of supersymmetric particles can be computed within DarkSUSY or with the help of 

5_( ' external programs such as FeynHiggs, ISASUGRA and SUSPECT. Accelerator bounds can 

C^ , be checked to identify viable dark matter candidates. DarkSUSY also computes a large 

variety of astrophysical signals from neutralino dark matter, such as direct detection in 

low-background counting experiments and indirect detection through antiprotons, an- 

tideuterons, gamma-rays and positrons from the Galactic halo or high-energy neutrinos 

from the center of the Earth or of the Sun. Here we describe the physics behind the 

package. A detailed manual will be provided with the computer package. 
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1 Introduction 

During the past few years, remarkable progress has been made in cosmology, both observa- 
tionally and theoretically. One of the outcomes of these rapid developments is the increased 
confidence that most of the mass of the observable Universe is of an unusual form, i.e., not 
made up of ordinary baryonic matter. Recent analyses combining high-redshift supernova 
luminosity distances, microwave background fluctuations and the dynamics and baryon frac- 
tion of galaxy clusters indicate that the present mass density of matter in the Universe ^m = 
pM /Pciit normahzed to the critical density pcrit — 3iJg/(87rG'Ar) = h? x 1.9 • 10~^^ gcm"'^ 
is 0.1 < ^Mh^ < 0.2, which is considerable higher than the value VtBh^ < 0.023 allowed by 
big bang nucleosynthesis ^. Here h ~ 0.7± 0.15 is the present value of the Hubble constant 
in units of 100 km s~^ Mpc~^. A recent addition to the wealth of experimental data which 
support the existence of non-baryonic dark matter is the WMAP microwave background re- 
sults Pj. In a joint analysis of the WMAP data together with other CMBR experiments, 
large-scale structure data, supernova data and the HST Key Project, the WMAP team report 
fleh^ = 0.0224 ± 0.009, nnh^ = 0.135 ± 0.009, and 17a = 0.73 ± 0.04. Subtracting the 
baryonic contribution flBh^ from the matter density VIm^? leaves a non-baryonic cold dark 
matter density QcDAih'^ = 0.113 ± 0.009. 

Also from the point of view of structure formation, non-baryonic dark matter seems to 
be necessary, and the main part of it should consist of particles that were non-relativistic at 
the time when structure formed (cold dark matter, CDM), thus excluding light neutrinos. 
Under reasonable assumptions, the WMAP collaboration, using also galaxy survey and Ly-a 
forest data, limit the contribution of neutrinos to n,yh^ < 0.0076 (95 % c.L). 

A well-motivated particle physics candidate which has the required properties is the light- 
est supersymmetric particle, assumed to be a neutralino O^E]. (For thorough reviews of 
supersymmetric dark matter, see [H1I3-) Although supersymmetry is generally accepted as a 
very promising enlargement of the Standard Model of particle physics (for instance it would 
solve the so-called hierarchy problem which essentially consists of understanding why the 
electroweak scale is protected against Planck-scale corrections), little is known about what 
a realistic supersymmetric model would look like in its details. Therefore, it is a general 
practice to use the simplest possible model, the minimal supersymmetric enlargement of the 
Standard Model (the MSSM), usually with some additional simplifying assumptions. Of 
course, there is no compelling reason why the actual model, if nature is supersymmetric at 
all, should be of this simplest kind. However, the MSSM serves as a useful template with 
which to test current ideas about detection, both in particle physics accelerators and in dark 
matter experiments and contains many features which are expected to be universal for any 
supersymmetric dark matter model. In fact, the knowledge gained by studying the MSSM in 
detail may be of even more general use, since it provides one specific example of a WIMP, a 
weakly interacting massive particle, which is generically what successful particle dark matter 
models require.^ 

Over several years, we have developed analytical and numerical tools for dealing with 
the sometimes quite complex calculations necessary to go from given input parameters in 
the MSSM to actual quantitative predictions of the relic density of the ncutralinos in the 
Universe, and the direct and indirect detection rates. The program package, which we have 
named DarkSUSY, has now reached a high level of sophistication and maturity, and we have 
released it publicly for the benefit of the scientific community working with problems related 
to dark matter. This paper describes the basic structure and the underlying physical and 
astrophysical formulas contained in DarkSUSY, as well as examples of its use. The version of 

^A completely different type of particle is the axion 8 , which is very light, but was never in thermal 
equilibrium. Its phenomenology is very different and will not be treated here. 



the package described in this paper is DarkSUSY 4.1. 

For download of the latest version of DarkSUSY and for a more technical manual, please 
visit the official DarkSUSY website, Ihttp : //www . p hysto . se/~eds j o/darksusy/| 

2 Definition of the Supersymmetric model 

We work in the framework of the minimal supersymmetric extension of the Standard Model 
defined by, besides the particle content and gauge couplings required by supersymmetry, the 
superpotential (the notation used is that of [H] which marked the beginning of the development 
of DarkSUSY, and is similar to \W\ ) 

W = £,, [-e*j,YEV^Hi - d*j,Y o^lHI + ^*RYuCi\Hi ~ ^^HlH^j (1) 

and the soft supersymmetry-breaking potential 

Vsott = ey (^-e*j,AEYEilHi - dl^AcYz^q^^iJ^' + u*j,AuYu^lH^ - B^HIH^ + h.c.) 
+Hi*mlHl + Hl*mlHi 

+ ^MiBB + ]-M2 {w^W^ + 2W+W-\ + ^r^hgg. (2) 

We give these and the following expressions since they contain our sign conventions. It should 
be noted that various authors use various sign conventions, and many errors, often difficult to 
find, can be avoided by keeping careful track of the signs, as we have tried to do consistently 
in DarkSUSY. Here i and j are SU(2) indices (ei2 = +1). The Yukawa couplings Y, the 
soft trilinear couplings A and the soft sfermion masses M are 3x3 matrices in generation 
space, e, 1, u, d and q are the superfields of the leptons and sleptons and of the quarks and 
squarks. A tilde indicates their respective scalar components. The L and R subscripts on 
the sfermion fields refer to the chirality of their fermionic superpartners. i?, W^ and W are 
the fermionic superpartners of the U(l) and SU(2) gauge fields and g is the gluino field. /^ is 
the Higgsino mass parameter. Mi, M2 and M3 are the gaugino mass parameters, B is a soft 
bilinear coupling, while m\ and m\ are Higgs mass parameters. 

These input parameters are contained in common blocks in the program. The full set of 
input parameters in version 4.1 of DarkSUSY, to be given at the weak scale, is m^, tan/3, 
/Z, Ml, M2, M3, AEaa. Auaa, Aoaa, M^^^, Ml^^, M^^^, M^^^, Af|„„ (with a = 1,2,3). 
The user may either provide these parameters directly to DarkSUSY or take advantage of 
the implementation of a MSSM pre-defined through a reduced number of parameters. In this 
model, the basic set of parameters is /x, M2, ttia, tan/3, mo. At, Af,. Here ttia is the mass 
of the CP-odd Higgs boson and tan/3 denotes the ratio, 'y2/'yi, of the vacuum expectation 
values of the two neutral components of the SU(2) Higgs doublets. The parameters mo, At 
and Ai, are defined through the simplifying Ansatz: Mg = M^/ = Mij = M^; = M^ = mol, 
Au = diag(0, 0, At), Ad = diag(0, 0, A), A^ = 0. 

Below we will give some details to clarify our convention and additional features. Rele- 
vant quantities for phenomenological studies, such as the particle masses and mixings, are 
consistently computed by DarkSUSY and available in arrays. The supersymmetry part of the 
program can thus be used for many applications, in particular for accelerator-based physics 
studies. Particle decay widths are also available, but currently only the widths of the Higgs 
bosons are calculated, the other particles having fictitious widths of 1 or 5 GcV (for the sole 
purpose of regularizing annihilation amplitudes close to poles) . 



2.1 Neutralino and chargino sectors 

The neutralinos x? are linear combinations of the superpartners of the neutral gauge bosons, 
B, W3 and of the neutral Higgsinos, Hi, i/^ ■ ^^ this basis, their mass matrix is given by 



-^1,2,3,4 






(3) 



where g and g' are the gauge coupling constants of SU(2) and U(l). S^^ and S44 are the 
most important one-loop corrections. These can change the neutralino masses by a few GeV 
up or down and are only important when there is a severe mass degeneracy of the lightest 
neutralinos and/or charginos. The expressions for S33 and S44 used in DarkSUSY are taken 
from ^2 EI (the tree-level values can optionally be chosen). 

The neutralino mass matrix, Eq. Q, is diagonalized analytically and evaluated numeri- 
cally to give the masses and compositions of four neutral Majorana states, 



r.O 



X 



= N,iB + N,2W^ + N^Hl + N,4H^^, (4) 



the lightest of which, Xj* to be called x for simplicity, is then the candidate for the particle 
making up the dark matter in the Universe. The neutralinos are ordered in mass such that 
TTi^o < 771^0 < TO^o < m^o and the eigenvalues are real with a complex N . 

The charginos are linear combinations of the charged gauge bosons W and of the charged 
Higgsinos H^ , H2 ■ Their mass terms are given by 

{W- Hr)M^± (Zl)+h.c. (5) 



H^ 



Their mass matrix, 



f M2 gv2 
\gvi ^ 



■A^x- = (r '?). (6) 



is diagonalized by the following linear combinations 

X7 = UaW- + U,2Hi, (7) 

X+ = VaW+ + V,2H+. (8) 

We choose det(U) = 1 and U*My± V^ = diag(TO -± , m-± ) with non-negative chargino masses 
"i^± > 0. We do not include any one-loop corrections to the chargino masses since they are 
negligible compared to the corrections ^aa and 644 introduced above for the neutralino masses 

nu. 

2.2 Sfermion masses and mixings 

When discussing the squark mass matrix including mixing, it is convenient to choose a basis 
where the squarks are rotated in the same way as the corresponding quarks in the Standard 
Model. We follow the conventions of the Particle Data Group ^^l and put the mixing in the 
left-handed d-quark fields, so that the definition of the Cabibbo-Kobayashi-Maskawa matrix 
is K = V1V2, where Vi (V2) rotates the interaction left-handed w-quark (d-quark) fields 



to mass eigenstates. For sleptons we choose an analogous basis, but since in DarkSUSY 4.1 
neutrinos are assumed to be massless, no analog of the CKM matrix appears. 
The general 6 x 6 u- and d-squark mass matrices are 

j^^^/Ml + raim^ + Dl^l mt ( A^, - /,* cot /3) \ 
" V (Ac/ - M cot /3)m„ M^ + m„mt + D^j^l / ' 

j^2^(^^MlK + maml + D'l^l ml{A^^ - ^^* t^n (3) \ 
"■ \ (Az)-/itan/3)md M], + m\md + D%j^l ) ' 

The general sneutrino and charged slepton mass matrices are (for massless neutrinos) 

Ml = Ml + Dl^l (11) 

j^2^(yil + xa,xnl + Dl^l mt(At, - ^* tan/3) \ ,^^. 

" \ (AB-^tan/3)me M| + mtnie + Dfj^jl / ' ^ ' 
Here 

d{l = ^l cos 2/3(T3/ - e/ sin^ Ow), (13) 

D^^ = m| cos(2/3)e/ sin^ % (14) 

where T3/ is the third component of the weak isospin and e/ is the charge in units of the 
absolute value of the electron charge, e. In the chosen basis, we have m„ = diag (tou, rn^, rrit), 
m.d = diag(md,ms,mb) and mg = diag(TOe,TO^,m^). 

The slepton and squark mass eigenstates fk {vk with k = 1,2,3 and Cfe, Uk and dk with 
k — 1, . . . , 6) diagonalize the previous mass matrices and are related to the current sfermion 
eigenstates /^^ and fjia (a = 1, 2, 3) via 

6 

ha = V/fe^>i^ (15) 



fe=l 



I Ra 



fe=l 



hn% (16) 



The squark and charged slepton mixing matrices Tul^ ^ur, Tdl, ^dr, ^el, and Ter have 
dimension 6x3, while the sneutrino mixing matrix F^^ has dimension 3x3. 

The current version of DarkSUSY allows only for diagonal matrices Ajj, Ad, Ae, Mq, 
M(7, Md, Mb, and M^. This ansatz, while not being the most general, implies the absence 
of tree-level flavor changing neutral currents in all sectors of the model. It also allows the 
squark mass matrices to be diagonalized analytically. For example, for the top squark one 
has, in terms of the top squark mixing angle 0^, 

rM = ^'dk = cos Oi, T^l = -r*}^ = sin e^. (i?) 

The sfermion masses are obtained with the diagonalization just described. 

To facilitate comparisons with the results of other authors, DarkSUSY allows for special 
values of the sfermion masses to be set in the program. A common value can be assigned to 
all squark masses, and an independent common value to all slepton masses. Alternatively, 
to enforce the sfermions to be heavier than the lightest neutralino (which we want to be 
the LSP), the squarks and sleptons masses can be set equal to the maximum between the 
neutralino mass and a specified value. In the special cases just described, no mixing is 
assumed between sfermions. It must be noted that the special choices described in this 
paragraph are mathematically inconsistent within the MSSM, but are often made for the 
sake of simplicity. 



2.3 Higgs sector and interface to FeynHiggs 

The Higgs masses receive radiative corrections and DarkSUSY includes several options for 
calculating these. The default in DarkSUSY is to use FeynHiggsFast J^ for the Higgs mass 
calculations. When higher accuracy is needed, it is possible to instead use the full FeynHiggs 
^iJ^ package. 

2.4 Interface to the mSUGRA codes ISASUGRA and SUSPECT 

Given the modular structure of DarkSUSY, the user may also run the package using as input 
for the MSSM definition the output, still at the low energy scale, from an external package. 
An example of usage under such mode is given in case of the minimal supergravity (mSUGRA) 
model: in the release, we provide an interface to the output of the ISASUGRA code, as included 
in ISAJET [T7I for the current ISASUGRA 7.69 version, as well as an an interface to SUSPECT 
[TH] which can be used as an alternative. The interfaces to DarkSUSY are at the level of the 
full spectrum of masses and mixings, including, for consistency, those for the Higgs sector. 
Of importance for relic density calculations near a Higgs boson resonance is the possibility 
of including supersymmetric corrections to the bottom, top, and tau Yukawa couplings as 
supplied by ISASUGRA. 

3 Experimental constraints 

Accelerator bounds can be checked by a call to a subroutine. By modifying an option, the 
user can impose bounds as of different moments in time. The default option in version 4.1 
adopts the 2002 limits by the Particle Data Group ^HI modified as described below. The 
user is also free to use his or her own routine to check for experimental bounds, in which case 
there is only need to provide an interface to DarkSUSY. 

For the theoretical prediction of the rare decay 6 — > 57 we have implemented the complete 
next-to-leading order (NLO) Standard Model calculation and the dominant NLO supersym- 
metric corrections. For the NLO QCD calculation of the Standard Model prediction we 
have used the expressions of reference j^Oj into which we have inserted the updated so-called 
"Magic numbers" of |2J. Our implementation of the Standard Model calculation gives a 
branching ratio BR[B ^ Xs 7] = 3.72 x 10"''^ for a photon energy greater than mb/20. 
This result agrees to within 1% with the result stated in [21], but is around 10% larger 
than the result of previous analyses. The latter is due to the fact that in the reference 
|2()) they replaced rnP°^'^ jmf^ ^ by m^ {^)/m^° ^ (with /i e [TTic^^b]) in the matrix element 
(Xs7 I (sc)y_A(c6)\/-A I &). 

The supersymmetric correction to 6 ^ 57 has been divided into a contribution from a 
two Higgs doublet model and a contribution from supersymmetric particles. The expressions 
for the NLO contributions in the two Higgs doublet model has been taken from ,22|. The 
corrections due to supersymmetric particles are calculated under the assumption of minimal 
flavour violation. The dominant LO contributions which are valid even in the large tan/3 
regime was taken from ref. [^Hl, and we also followed their guideline on how the NLO QCD 
expressions of [23 should be expanded to the large tan /3 regime. 

For the current experimental bound on & — > 57 we take the value stated by the Particle 
Data Group 2002 ^Hl- This is an average between the CLEO and the Belle measurements 
and amounts to BR[i? — > Xs 7] = (3.3 ±0.4) x lO^''. To this we add a theoretical uncertainty 
which we set to ±0.5 x 10^"*. The final constraint on the branching ratio then becomes 
2.0 X 10-"* < BR[B -^Xsi\< 4.6 x 10"''. 



Recently, much interest has been given to the possible contribution of supersymmetry 
to {g — 2)p. Although the discrepancy with the Standard Model result is now below 3a, 
we include for convenience a calculation of the anomalous moment of the niuon (g — 2)^ in 
DarkSUSY. 

4 Calculation of relic density 

The WMAP microwave background experiment [2| , combined with other sets of data, gives 
a quite precise determination of the cold dark matter density flcDAih'^ = 0.113 ± 0.009. 
We would like DarkSUSY to compute the relic density of neutralinos to at least the same 
precision. 

We use in DarkSUSY the full cross section and solve the Boltzmann equation numerically 
with the method given in |25l f^ . In this way we automatically take care of thresholds and 
resonances. 

When any other supersymmetric particles are close in mass to the lightest neutralino 
they will also be present at the time when the neutralino freezes out in the early Universe. 
When this happens so-called coannihilations can take place between all these supersymmetric 
particles present at freeze-out. Coannihilations were first pointed out by Binetruy, Girardi, 
and Salati [571 in a non-supersymmetric model with several Higgs bosons. Griest and Seckel 
PHj investigated them in the MSSM for the case where squarks are of about the same mass as 
the lightest neutralino. Later, coannihilations between the lightest neutralino and the lightest 
chargino were investigated in [211 I30L [TT] . Several authors have also included coannihilations 
with sfermionsSl, 32, 33, 3i EHl EE EZI [38 , 39, 40, 41, 42, 43, 44,!iS||iS|- In DarkSUSY 
we have implemented all coannihilations between the neutralinos, charginos and sfermions 
as calculated in 146;. 

Compared to other recent calculations, we believe ours is the most precise calculation 
available at present. The standard lore so far has been to calculate the thermal average of 
the annihilation cross section by expanding to first power in temperature over mass and im- 
plementing an approximate solution to the evolution equation which estimates the freeze out 
temperature without fully solving the equation (see, e.g., Kolb and Turner 0I|). Sometimes 
this is refined by including resonances and threshold corrections PSI. Among recent studies, 
this approach is taken in e.g. Refs. [3^1 1^. Other refinements include, e.g., solving the 
density evolution equation numerically but still using an approximation to thermal effects in 
the cross section [21 021 OSl 1^ OH] , or calculating the thermal average precisely but using 
an approximate solution to the density equation :38^. >39l BUj . At the same time, only in a few 
of the quoted papers the full set of initial states has been included. As already mentioned, 
the present calculation includes all initial states, performs an accurate thermal average and 
gives a very precise solution to the evolution equation. Though the inclusion of initial state 
sfermions in the DarkSUSY package is a new feature introduced in the present work, other 
groups J41ll42ll^ have earlier introduced some sfermion coannihilations in an interface with 
the old DarkSUSY version. 

To gain calculational speed we only include the particles with masses below fco^n^. The 
mass fraction parameter fco is by default set to 2.1 or 1.4 depending on how the relic density 
routines are called (very high precision or fast calculation), but can be set to any value by 
the user. 

4.1 The Boltzmann equation 

We will here outline the procedure developed in [5^1 which is used in DarkSUSY. For more 
details, see |2E1- 



Consider annihilation of N supersymnietric particles with masses rm and internal degrees 
of freedom gi. Order them such that mi < ni2 < • • • < ttt-n-i < rriN. For the lightest 
neutralino, the notation mi and m-^ will be used interchangeably. 

Since we assume that i?-parity holds, all supersymmetric particles will eventually decay 
to the LSP and we thus only have to consider the total number density of supersymmetric 
particles n = X]i=i"'i- Furthermore, the scattering rate of supersymmetric particles off 
particles in the thermal background is much faster than their annihilation rate, because 
the scattering cross sections ct^j,- are of the same order of magnitude as the annihilation 
cross sections Cy but the background particle density nx is much larger than each of the 
supersymmetric particle densities Ui when the former are relativistic and the latter are non- 
relativistic, and so suppressed by a Boltzmann factor 0S]. Hence, the Xi distributions remain 
in kinetic thermal equilibrium during their freeze-out. Combining these effects, we arrive 
at the following Boltzmann equation for the summed number density of supersymmetric 
particles 

— = -3Hn - {a,sv) (n^ - n^^) (18) 

where 

KW=^Kt;„)^^. (19) 



u 



with 



iPi-PjY 



Ziyy^J, 



m,7m 



3 



E^,E, 



(20) 



4.2 Thermal averaging 



Using the Maxwell-Boltzmann approximation for the velocity distributions one can derive 
the following expression for the thermally averaged annihilation cross section |26| 

/o°° dp,splaW,«Ki (^) 

{<^effV) = ^ . (21) 

mJT E,— ^^2(^) 

where Ki {K2) is the modified Bessel function of the second kind of order 1 (2), T is the 
temperature, s is one of the Mandelstam variables and 






In this equation, we have defined the momentum 



[s - (m, + mjf] ^'"^ [s - (m, - m.jf\ ^^^ 
P. = Vi ' ^ ^ 



the invariant annihilation rate 



w^.. 



^PijVscTij = 4:aijJ{pi ■ pjY - mjm'j ^ 4E.iEjaijVij (24) 



and the effective momentum 



f 



Pcff = Pu 



imi 



(25) 



Since Wij{s) ~ for s < {nii + rrij)'^, the radicand in Eq. H22|) is never negative. For a 
two-body final state, Wij is given by 



W, 



2-body 



»J 



l6TT^gtgjSf^/s 



E 



internal d.o.f. 



\Mfdn, 



(26) 



where k is the final center-of-mass momentum, Sf is a symmetry factor equal to 2 for identical 
final particles, and the integration is over all possible outgoing directions of one of the final 
particles. As usual, an average over initial internal degrees of freedom is performed. 

4.3 Annihilation cross sections 

In DarkSUSY, all two-body final state cross sections at tree level are computed for all coan- 
nihilations between neutralinos, charginos and sfermions. A complete list is given in table 13 
in Appendix IXI 

The calculation of the relic density is, due to its complexity, the most time-consuming 
task of DarkSUSY. For the neutralino-neutralino, chargino-neutralino and chargino-chargino 
initial states, to achieve efficient numerical computation, contributing diagrams have been 
classified according to their topology (s-, t- or w-channel) and to the spin of the particles 
involved. The helicity amplitudes for each type of diagram have been computed analytically 
with Reduce |49j using general expressions for the vertex couplings. The Reduce output has 
been automatically converted to FORTRAN subroutines called by DarkSUSY. 

The strength of the helicity amplitude method is that the analytical calculation of a given 
diagram only has to be performed once and the summing of the contributing diagrams for 
each given set of initial and final states can be done numerically afterwards. 

The Feynman diagrams are summed numerically for each annihilation channel ij —^ kl. 
We then sum the squares of the helicity amplitudes and sum the contributions of all annihi- 
lation channels. Explicitly, DarkSUSY computes 



dWcs _ s-y PijPki sr-y 
dcose "^32^5fe,yi^ ,^. 

tjKL nclicitics 



Y^ M{ij ^ kl) 



diagrams 



(27) 



where 9 is the angle between particles k and i. Finally, a numerical integration over cos 9 is 
performed by means of an adaptive method (SOI . 

In rare cases, there can be resonances in the t- or w-channels. This can happen when 
the masses of the initial state particles lie between the masses of the final state particles. 
At certain values of cos0, the momentum transfer is time-like and matches the mass of the 
exchanged particle. We regulate the divergence by assigning a small width of 5 GeV to the 
neutralinos and charginos and 1 GeV to the sfermions. The results arc not sensitive to the 
exact choice of this width. 

For the coannihilation diagrams with sfermions, the calculations are done with Form |51| 
and automatically converted into FORTRAN subroutines. 

In the relic density routines, the calculation of the effective invariant rate Weff is the most 
time-consuming part. Fortunately, thanks to the remarkable feature of Eq. H21|l . Weff(Poff) 
does not depend on the temperature T, and it can be tabulated once for each model. 

To perform the thermal average in Eq. H21(l , we integrate over p^s by means of adaptive 
gaussian integration, using a spline routine to interpolate in the {peS,Wes) table. To avoid 
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numerical problems in the integration routine or in the spline routine, we split the integration 
interval at each sharp threshold. We also explicitly check for each MSSM model that the 
spline routine behaves well close to thresholds and resonances. 

Wc then integrate the density evolution equation, Eq. H18|) . For numerical reasons, we 
do not integrate the equation directly, but instead rephrase it as an evolution equation for 
the abundance, Y = n/s (with s being the entropy density) and use x = m^/T as our 
integration variable instead of time (see |26j for details) . The numerical integration is subtle 
since the equation is "stiff." For this purpose, we developed an implicit trapezoidal method 
with adaptive stepsize. In short, if the equation for Y is written as dY/dx = X{Y'^ — Y^), 
with Ye being the equilibrium value at temperature T, the numerical integration is based on 
the recurrence 

Yt+1 = , ' (28) 

1 + VI + hx,+ia 

where h is the stepsize and Ci = 2Yi + h[Xi-^-iY^^^-^ + Xi{Y^^^ — Y-^)]. The stepsize is 
adapted (reduced) if \{Yi^i — yi+i)/Yi^i\ exceeds a given tolerance. Here y^+i = (c,;/2)(l + 
^/TThX^^)-^ with c, = 4(r, + hX^Y^,). 

There are some loop-induced final states, such as two gluons, two photons or a Z" boson 
and a photon which could in principle give a non-negligible contribution to the annihilation 
rate and lower the relic abundance somewhat. The cross sections for these processes are 
very complicated already in the limit of zero velocity, and we therefore assume that the 
invariant rate W for these one-loop processes are constant and equal to their zero-momentum 
expressions. These processes can be excluded from the calculation by setting appropriate 
parameters in the code. 

The relic density routines can be called in a precise way where all integrations are per- 
formed to a precision better than 1% and coannihilations are included up to a mass difference 
of fco = 2.1. It can also be called in a faster way, where the precision of the integrations are 
of the order of 1% and coannihilations are included up to fco — 1-4. Usually, the difference 
between the precise and fast method is completely negligible, but in rare cases it can be of a 
few %. The fast option is considerably faster and should be sufficient for most purposes. For 
advanced users, it is also possible to manually decide exactly which coannihilation processes 
to include. 

For users with less demand of calculational precision, we also provide in DarkSUSY the 
option of a "quick-and-dirty" method of computing the relic density, essentially according to 
the textbook treatment in [Uj. It should be realised, however, that this may sometimes give 
a computed relic density which is wrong by orders of magnitude. 

4.4 A note about degrees of freedom and the annihilation rate 

We end this section with a comment on the internal degrees of freedom gi. A neutralino 
is a Majorana fermion (it is its own antiparticle) and has two internal degrees of freedom, 
g^o =2. A chargino can be treated either as two separate species xt ^^'^ Xi~i each with 
internal degrees of freedom g^+ — g^- — 2, or, more simply, as a single species Xi with 
g ± = 4 internal degrees of freedom. We choose the latter convention in DarkSUSY and use 

the analogous conventions for sfermions (see j33 for details). 

The counting of states in annihilation processes for Majorana fermions is non-trivial, 
and has led to a factor-of-two ambiguity in the literature which also propagated into earlier 
versions of DarkSUSY. The current version of DarkSUSY contains the correct normalization 
of annihilation rates, namely av in the Boltzmann equation and av/2 for annihilation in the 
halo. The clearest way to see the origin of the factor of 1/2 is probably to go back to the 
Boltzmann equation, |52| . In essence, one can view av as the thermal average (averaged over 
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momentum and angles) of the cross section times velocity in the zero momentum limit; in this 
average one integrates over all possible angles. For identical particles in the initial state, one 
includes each possible initial state twice, therefore one needs to compensate by dividing by a 
factor of 2; the prefactor in the zero-momentum limit becomes then av/2. In the Boltzmann 
equation describing the time evolution of the neutralino number density the 1/2 does not 
appear as it is compensated by the factor of 2 one has to include because 2 neutralinos are 
depleted per annihilation, but we need to include the factor of 1/2 explicitly in other cases 
where we need the annihilation rate (like for annihilation in the halo). Another way to view 
this is to think of <t as to the annihilation cross section for a given pair of particles. Let 
the number of neutralinos in a given volume be N] the annihilation rate would be given by 
av times the number of pairs, which is N{N — l)/2. In the continuum limit this reduces to 
avn? /2. 

5 Halo models 

The modelling of the distribution of dark matter particles in the Milky Way dark halo plays 
a major role in estimates of dark matter detection rates. On this issue, however, there is 
no well-established framework we can refer to. Available dynamical measurements, such 
as, e.g., the mapping of the rotation curve, the local field of acceleration of stars in the 
direction perpendicular to the disk, or the motion of the satellites in the outskirts of the 
Galaxy, provide some constraints on the dark matter density profile, but lack the precision 
needed to derive a refined model. The guideline for the future will be N-body simulations of 
hierarchical clustering in cold dark matter cosmologies, which are starting to resolve the inner 
structure of individual galactic halos. At present, however, the translation of such results 
into the detailed model we need for the Milky Way still relies on large extrapolations, and, 
to some extent, faces the problem of possible discrepancies between some of its properties 
and observations in real galaxies (for a recent review see, e.g., j53p. 

In light of these large uncertainties, the definition of the model for the dark matter halo 
in DarkSUSY has been kept in a very general format. Two sets of properties need to be 
implemented: 

a) Local properties: To derive counting rates in direct detection and x-induced neutrino 
fluxes from the center of the Sun and the Earth, the user has to specify the halo 
mass density po at our galactocentric distance Rq and the relative particle velocity 
distribution /(m), where u is the modulus of the velocity u in the rest frame of the Sun 
(or Earth, respectively) and an average of incident angles has been assumed. Options 
to set f{u) include the possibility to implement the expression valid for a (truncated) 
isothermal sphere, or the numerical interpolation of a table of pairs {fi,Ui) provided 
by the user. We have also implemented routines which compute the function /(u) in 
the Earth or the Sun rest frame for a given isotropic velocity distribution function in 
the Galactic reference frame as needed for direct detection and for neutralino capture 
rates by the Sun. Finally a numerical table for f{u) in the Earth frame and including 
the modelling of the diffusion process of neutralinos through the solar system |^ is 
available to compute capture rates by the Earth. 

b) Global properties: To compute indirect signals from pair annihilations in the halo, the 
full dark matter mass density profile is needed. Charged cosmic ray fluxes are com- 
puted in two dimensional models and a generic axisymmetric profile p{R, z) can be 
implemented by the user. Line of sight integrals with angular averaging over accep- 
tance angles, needed to compute gamma-ray fluxes, are given for spherically symmetric 
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profiles p{r). The options to specify p include several analytic forms, with most classes 
of profiles proposed in recent years |55ll56[IK7ll58| . as well as the numerical interpola- 
tion of a table of pairs {pi, Vi) provided by the user. It is also possible to switch on an 
option to compute the signals for annihilations taking place in a population of small, 
unresolved clumps. In this case the user should specify the probability distribution for 
the clumps and the average enhancement of the source per unit volume compared to 
the smooth halo scenario. 

An option to rescale po a-nd p{r) by the quantity Q,h^ j Q.^i^h? , where Q.h? is the neutralino 
relic abundance and ilmin/i^ a minimum reference value for neutralinos to provide most of 
the dark matter in our Galaxy, is available for the case where subdominant dark matter 
candidates are considered. The effect can be introduced a posteriori as well for all detection 
techniques except for neutrino fluxes, where the response to rescaling is non-linear. 

A separate package interfaced to DarkSUSY with self consistent pairs {p{r),f{u)) in 
ACDM inspired models which fit available dynamical constraints will be available shortly 
from one of the authors |59| . 

6 Detection rates 

For the detection rates of neutralino dark matter we have used the rates as calculated in Rcfs. 
|52l 1601 El Ifill lf)2l 1681 1641 1651 ES| , with some improvements which we report in this Section 
where the formulas used in DarkSUSY are presented. 

6.1 Direct detection 

The current version of DarkSUSY provides the neutralino-proton and neutralino-neutron scat- 
tering cross sections (spin-independent and spin-dependent), as well as cross sections and form 
factors for elastic scattering of neutralinos off nuclei. 

The rate for direct detection of galactic neutralinos can be written as 

^=Vc,^^^^i^#l^/ /(M,3,. (29) 

dE ^ 2m^Mix A>^m,b/2mJ. -v 

The sum is over the nuclear species in the detector, ci being the detector mass fraction 
in nuclear species i. Mi is the nuclear mass, and p^i — m^Mi/{m^ + Mi) is the reduced 
neutralino-nucleus mass. Moreover, p^ is the local neutralino density, v the neutralino 
velocity relative to the detector, v = |v|, and f{w,t)d?v the neutralino velocity distribution. 
Finally, a^i is the total scattering cross section of a WIMP off a fictitious point-like nucleus, 
\Fi{qi)\'^ is a nuclear form factor that depends on the momentum transfer qi = \/2MiE and 
is normalized to Fi(0) = 1. The integration is over all neutralino speeds that can impart a 
recoil energy E to the nucleus. 

The cross section a^i scales differently for spin-dependent and spin-independent interac- 
tions. For spin-independent interactions with a nucleus with Z protons and A — Z neutrons, 
one has 

all = ^\ZGP + (A -Z)G':\\ (30) 

where Gf and G" are the scalar four-fermion couplings of the WIMP with point-like protons 
and neutrons, respectively (see, e.g., |S7], and below). As default, the spin-independent form 
factor in DarkSUSY is taken to be of the Helm form 
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with ji a spherical Bessel function of first kind, Ri = y/W~~5^, R = [0.9{M/GeV)^/^ + 
0.3] fm and s ~ 1 fm. An exponential form factor is also available as an option. 
For spin-dependent interactions, one has instead, at zero momentum transfer, 



SD ^^^Y^ J + ^ 



SD ^ j:^^_Us,)Gp + {Sr.)G:\\ (32) 



where J is the nuclear spin, (Sp) and {S„) are the expectation values of the spin of the 
protons and neutrons in the nucleus, respectively, and G^ and GJJ are the axial four-fermion 
couplings of the WIMP with point-like protons and neutrons |()7I l()8| . At finite momentum 
transfer, the spin-dependent cross section times the form factor reads 

4u2. 
^^f l^'°('Z)P = ^jfY m?Spp{q) + {G-fSnn{q) + GIG-Sp^iq)] , (33) 

where Spp{q) = 5*00(9) + Sii{q) + S'oi(q), Snn{q) = <5'oo(<7) + S'ii(<7) - Sfii{q), and Spn{q) = 
2[S'oo('z) — Sii{q)]. The spin structure functions 500(9), 5ii(g), and 5oi(<z) are defined in [5^ 
and are given in the literature. 

For protons and neutrons, the previous expressions reduce to 

crl = ^\Gl\\ af^ = ^\G:\\ a^° = ^|GSP, a^;^ ^ ^-^\GZ? . (34) 

For the neutralino, the scalar and axial four-fermion couplings with the proton and neu- 
tron arise from squark, Higgs and Z boson exchange. In DarkSUSY, the default expressions 
for a nucleon N = n,p are 



■)' 



(6 2 2 \ 

gzxx.9^.. ^ 1 y 9l^.x, +9R.1.X, ^ (35) 

8 ^1 m^, J 

where gabc are the coupling constants in the vertex involving particles abc (see |3] and |65| 
for explicit expressions). The more complicated expressions of Drees and Nojiri [TOj are also 
available as an option. 

Default values of the parameters in DarkSUSY are [^[72 (with {N\qq\N) = Jto'^^n /rriq) 

/P„ = 0.023, /?, = 0.034, /l^, = 0.14, /|^^ = /|:, = /^^ = 0.0595, (37) 

/?„ = 0.019, /?, = 0.041, /^, -0.14, /^, = /^, - /^', = 0.0592. (38) 

(Au)p = (Ad)„ = 0.77, (Ad)p = (Au)„ = -0.40, (As)p = (As)„ - -0.12. (39) 

Other sets of values for these parameters are available. These values can also be overridden 
by the user. 

6.2 Indirect detection 

There are several ways of detecting dark matter particles indirectly. Pair annihilations of dark 
matter particles in the Galactic halo may give an exotic component in positron, antiproton or 
antideuteron cosmic-rays and gamma-rays. There may also be annihilation in astrophysical 
environments where the density of neutralinos may be enhanced, such as annihilation in the 
center of the Earth or Sun (detected in neutrino telescopes through a high-energy neutrino 
fiux) or near the central black hole of the Galaxy. 
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6.2.1 Monte Carlo simulations 

In several of the indirect detection processes below we need to evaluate the yield of different 
particles per neutralino annihilation. The hadronization and/or decay of the annihilation 
products are simulated with Pythia "731 6.154 in essentially the same way (with a few dif- 
ferences) for all these processes and we here describe how the simulations are done. We can 
divide the simulations into two groups: a) annihilation in the Earth and the Sun and b) 
annihilation in the halo. In both cases the simulations are done for a set of 18 neutralino 
masses, m^ = 10, 25, 50, 80.3, 91.2, 100, 150, 176, 200, 250, 350, 500, 750, 1000, 1500, 2000, 
3000 and 5000 GeV. We tabulate the yields and then interpolate these tables in DarkSUSY. 
These two groups of simulations differ slightly in other aspects, namely 

a) Annihilation in the Earth and the Sun. In this case we are mainly interested in the flux 
of high energy muon neutrinos and neutrino-induced muons at a neutrino telescope. 
We simulate 6 'fundamental' annihilation channels, cc, 66, ti, t^t^ , W'^W~ and Z'^Z^ 
(if kincmatically allowed) for each mass listed above. The lighter Icptons and quarks do 
not contribute significantly and can safely be neglected. Pions and kaons get stopped 
before they decay and are thus made stable in the Pythia simulations so that they do 
not produce any neutrinos. For annihilation channels containing Higgs bosons, we can 
calculate the yield from these fundamental channels by letting the Higgs bosons decay 
in flight (see below). We also take into account the energy losses of B-mesons in the 
Sun and the Earth by following the approximate treatment of 'J4} but with updated 
_B-meson interaction cross sections as given in |65) . Neutrino- interactions on the way 
out of the Sun are simulated with Pythia including neutral current interactions and 
charged-current interactions as a neutrino- loss. The neutrino- nuclcon charged current 
interactions close to the detector are also simulated with Pythia and finally the mul- 
tiple Coulomb scattering of the muon on its way to the detector is calculated using 
distributions from ^^j ■ For more details on these simulations see [ZHl HE] ■ 

For each annihilation channel and mass we simulate 1.25 x 10^ annihilations and tabu- 
late the final results as a neutrino-yield, neutrino-to-muon conversion rate and a muon 
yield differential in energy and angle from the center of the Sun/Earth. We also tab- 
ulate the integrated yield above a given threshold and below an open-angle 9. We as- 
sumed throughout that the surrounding medium is water with a density of 1.0 g/cm"^. 
Hence, the neutrino-to-muon conversion rates have to be multiplied by the density of 
the medium. In the muon fluxes, the density cancels out (to within a few percent). For 
the neutrino- nucleon cross sections, we have used the parameterizations in |65j . 

b) Annihilation in the halo. The simulations are here simpler since we do not have a 
surrounding medium that can stop the annihilation products. We here simulate for 
8 'fundamental' annihilation channels cc, bb, ti, t^t~, W~^W~, Z'^Z'^, gg and /i^/i^. 
Compared to the simulations in the Earth and the Sun, we now let pions and kaons 
decay and we also let antineutrons decay to antiprotons. For each mass we simulate 
2.5 X 10^ annihilations and tabulate the yield of antiprotons, positrons, gamma rays 
(not the gamma lines), muon neutrinos and neutrino-to-muon conversion rates and 
the neutrino-induced muon yield, where in the last two cases the neutrino-nucleon 
interactions has been simulated with Pythia as outlined above. 

With these simulations, we can calculate the yield for any of the above mentioned particles 
for a given MSSM model. For the Higgs bosons, which decay in flight, an integration over 
the angle of the decay products with respect to the direction of the Higgs boson is performed. 
Given the branching ratios for different annihilation channels it is then straightforward to 
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compute the muon flux above any given energy threshold and within any angular region 
around the Sun or the center of the Earth. 

6.3 Neutrinos from the Sun and Earth 

One of the most promising indirect detection method jZJ] relies on the fact that scattering 
of halo neutralinos by the Sun and the planets, in particular the Earth, during the several 
billion years that the Solar system has existed, will have trapped these neutralinos within 
these astrophysical bodies. Being trapped within the Solar or terrestrial material, they 
will sink towards the center, where a considerable enrichment and corresponding increase of 
annihilation rate will occur. 

Searches for neutralino annihilation into neutrinos will be subject to extensive experi- 
mental investigations in view of the new neutrino telescopes (AMANDA, IceCube, Baikal, 
NESTOR, ANTARES) planned or under construction [75]. A high-energy neutrino signal in 
the direction of the center of the Sun or Earth is an excellent experimental signature which 
may stand up against the background of neutrinos generated by cosmic-ray interactions in 
the Earth's atmosphere. 

The detector energy threshold for "small" neutrino telescopes like Baksan, MACRO and 
Super-Kamiokande is around 1 GeV. Large area neutrino telescopes in the ocean or in Antarc- 
tic ice typically have thresholds of the order of tens of GeV, which makes them sensitive 
mainly to heavy neutralinos (above 100 GeV) |Z21- In DarkSUSY, the low energy cut-off can 
be set. 

Final states which give a hard neutrino spectrum (such as heavy quarks, r leptons and 
W OT Z bosons) are usually more important than the soft spectrum arising from light quarks 
and gluons. 

Neutralinos are steadily being trapped in the Sun or Earth by scattering, whereas anni- 
hilations take them away. The balance between capture and annihilation has the solution for 
the annihilation rate implemented in DarkSUSY 

Cc i__, 2 / t 



r^^_LtanhM-), (40) 



where the equilibration time scale r — l/\/CcCa, with Cc being the capture rate and Ca 
being related to the annihilation efficiency. In most cases for the Sun, r is much smaller 
than a few billion years, and therefore equilibrium is often a good approximation {N{t) = 0). 
This means that it is the capture rate which is the important quantity that determines the 
neutrino flux. For the Earth, r is, on the other hand, usually of the same order as, or much 
larger than, the age of the solar system, and equilibrium has often not occurred. In either 
case, in the program we keep the exact formula (|40|l (except for the modifications needed for 
a Damour-Krauss population of WIMPs, discussed below). 

For the actual capture rate calculations we have several expressions implemented in Dark- 
SUSY. As a default, we use the full expressions given in appendix A of |8(J| where we numer- 
ically integrate over the velocity distribution. To speed-up the calculations, it is possible to 
perform this integration only once and use a saved tabulated version for subsequent calls. 
In the capture rate calculations we also need the density profile of the Earth and the Sun 
and the chemical composition as a function of radius. For the Sun we use the solar model 
BP2000 |5T] , complemented with the estimates of the mass fractions of the heavier elements 
from |S^. For the Earth, we use the density profile of 83 with the compositions given in 
^^ (see |S3] for a table of these). For comparison, the approximate capture rate expressions 
in |H] are also available. 
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The capture rate induced by scalar (spin-independent) interactions between the neutrah- 
nos and the nuclei in the interior of the Earth or Sun is the most difficult one to compute, since 
it depends sensitively on the Higgs mass, form factors, and other poorly known quantities. 
However, this spin-independent capture rate calculation is the same as for direct detection 
treated in Section fo.ll Therefore, there is a strong correlation between the neutrino flux ex- 
pected from the Earth (which is mainly composed of spin-less nuclei) and the signal predicted 
in direct detection experiments 79, 85 . 

Due to the low counting rates for the spin-dependent interactions in terrestrial detectors, 
high-energy neutrinos from the Sun constitute a competitive and complementary neutralino 
dark matter search. Of course, even if a neutralino is found through direct detection, it will 
be extremely important to confirm its identity and investigate its properties through indirect 
detection. In particular, the mass can be determined with reasonable accuracy by looking at 
the angular distribution of the detected muons 86, 87 . 

The capture rate in the Earth is dominated by scalar interactions, where there may 
be kinematic and other enhancements, in particular if the mass of the neutralino almost 
matches one of the heavy elements in the Earth. As shown by Gould (SS], the Earth does 
not dominantly capture WIMPs from the halo directly. Instead, the Earth captures WIMPs 
that, due to gravitational interactions, have diffused around and become bound to the solar 
system. However, solar depletion of these bound WIMPs could be an important effect )89) . 
and as a default, DarkSUSY uses a new estimate of the velocity distribution in the solar 
system, where these solar depletion effects have been included jSJ. It is possible to change 
to a standard Gaussian distribution if the user prefers. 

There is also a possibility that there exists a special population of WIMPs, the Damour- 
Krauss population [^01, arising from WIMPs that have just skimmed the Sun's surface. As 
an optional choice, this population can be included in the calculation SQf . The enhancement 
caused by the new population is only important for a neutralino mass less than 150 - 170 
GeV (the exact number depending on details about the angular momentum distribution 1661). 
The total capture rate is computed according to the formulas in 66 , which take into account 
that the annihilation rates from the Earth will in general depend on time in a different way 
than the simple result in Eq. I|40|l . 

6.4 Indirect searches through antimatter signals 

Pair annihilation of ncutralinos in the Galactic halo produces the same amounts of matter and 
antimatter. As antimatter seems to be scarce in the Universe, apparently with no standard 
primary sources, there is a chance that by measuring antimatter in cosmic ray fluxes one 
may find evidence of the existence of dark matter particles (see, e.g., |^ |521 ESI US 1^ 
I96| \97[ l??HI IH9| ). In the current release of the code we consider neutralino induced fluxes 
of antiprotons, positrons and antideuterons. To produce estimates of such fluxes, there are 
several steps one needs to follow: i) Evaluate for each dark matter candidate the probability 
for a pair annihilation to take place locally in space, i.e. compute 1/2 (p^(x)/m^)'^aiinnV',^ 
a) Estimate the production rate of a given species by folding together, for each model, the 
branching ratios for the annihilation into a given two-body final state with the Monte Carlo 
simulation of the hadronization and/or decay of that state, as described in Section 16.2.11 
(except for D sources where we have implemented the prescription suggested in Ref. ,98, to 
convert from the p yield); Hi) Propagate these sources through the Galactic magnetic fields 
to make predictions for the induced cosmic ray fiuxes at the Sun's location in the Galaxy; iv) 
Include the effect of solar modulation to propagate the fiuxes from the interstellar medium 
to our location inside the solar system. The implementation in DarkSUSY is written in a 

*The origin of the factor of 1/2, missing in earlier versions of DarkSUSY, is explained in Sect. l4!4l 
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modular format which aUows the user to eventually modify and/or replace each of the four 
steps above independently. 

At tree level, the relevant final states of neutralino pair annihilations for p and D produc- 
tion are qq, W+W' , Z°Z°, W+ R- , ZH^, ZH^, H^H^ and H^Hl; among the qq states we 
have included in DarkSUSY all the heavier quarks (c, b and t). In addition, we have included 
the Z7 (ii'OO'') and the 2 gluon f ll01| : |102| ') final states which occur at one loop-level. The 
same list of final states is implemented for positrons, with the addition of i'^l^ . The e+e~ 
final state gives rise to a positron monochromatic source (however this is negligible in most 
particle physics setups, and tends to be smeared out by propagation). 

To model the propagation we have considered a semi-empirical diffusion equation |1U3[ 
I1U4| . in the steady state approximation, applied to a two-dimensional model with cylindrical 
symmetry and with free escape boundary conditions. The parameters in such a model should 
be fixed in agreement with values inferred from available cosmic ray data in analogous propa- 
gation models (see, e.g., the Galprop package |l()5j V For simplicity, rather than considering 
the most general setup, we restrain to cases in which, still including all relevant effects, the 
Green function of the transport equation can be derived analytically, so that we can avoid 
the CPU time-consuming problem of having to solve a partial differential equation for each 
dark matter candidate. Therefore, we do not include reacceleration effects but mimic them 
through a diffusion coefhcient which has a broken power law in rigidity as functional form. 
Also , for antiprotons jHE] and antideuterons we neglect energy losses (whenever a scattering 
with a nucleus takes place the particle is removed), while for positrons (HZ] we consider an 
average over space of the energy loss effect due to inverse Compton on starlight and the 
cosmic microwave background, and in addition the synchrotron radiation from the effects of 
the galactic magnetic field. For comparison, we allow also the option to treat the propagation 
of antiprotons according to the propagation models by Chardonnet et al. [H^ ^^d Bottino ct 
al. IHS], while for the positron we have implemented the option to use the leaky-box treat- 
ment given by Kamionkowski and Turner |93j or the numerical Green functions derived by 
Moskalenko and Strong 106 with the Galprop code |105j (the latter two cases however 
cannot be interfaced to a generic axisymmetric dark matter density profile, as for our default 
propagation model). 

Given a set of parameters for the propagation model, and a given neutralino number 
density profile or a given probability distribution of small clumps, the effect of propagation, in 
the cases we have considered, can be factorized out into effective energy-dependent "diffusion 
times", Tp(T), tq(T) and Tf^+{T), which are independent of the parameters defining the 
particle physics setup for the dark matter candidate. In some cases, such as when considering 
very large samples of neutralino candidates or when implementing singular halo profiles 
for which the computation of the diffusion times can be very CPU time-consuming, it is 
advantageous to exploit the option provided by DarkSUSY to tabulate the diffusion times 
over a given energy range (optionally saving the tabulation to disk for later use) and use an 
interpolation between tabulated values, rather than linking directly to the computation for 
each dark matter candidate. 

Finally, regarding solar modulation, we implement the one parameter model based on the 
analytical force-field approximation by Gleeson & Axford |1U7| for a spherically symmetric 
model. This approach is expected to be shghtly less accurate than, e.g., the full numerical 
solution of the propagation equation of the spherically symmetric model (^IHl), but again it 
is much less CPU time-consuming. DarkSUSY allows an output with both solar modulated 
and local interstellar fluxes, and the latter can eventually be solar modulated by the user 
with more sophisticated methods. For the positrons we allow for the option to reduce the 
effects of solar modulation by considering the positron fraction, e+/(e+ + e^), rather than the 
absolute positron fluxes. In this case an estimate of the background is needed: in DarkSUSY 
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we provide background e+ and e fluxes taken from |106| . 

6.5 Gamma rays 

Gamma-rays have a low enough cross section on gas and dust that the Galaxy is essentially 
transparent to them (except perhaps in the innermost part, very close to the region where 
a massive black hole is inferred); absorption by starlight and infrared background becomes 
efficient only for very far away sources and high-energy photons. 

The bulk of the gamma-rays from neutralino annihilations arise in the decay of neutral pi- 
ous produced in the fragmentation processes initiated by tree level final states |1()9I EOl IllOj 
(analogously to the other halo signals, in DarkSUSY we include all tree level final states 
and make use of a Monte Carlo simulation for fragmentation and decay processes, see Sec- 
tion |^^3J- However, vr*' production is common also to other astrophysical processes, and 
this may turn out to be a limiting factor to disentangle dark matter sources. At the same 
time, however, a relevant gamma-ray contribution may arise directly (at one-loop level) in 
two body final states; although such photons are much fewer than those from tt" decays they 
have a much better signature: neutralinos annihilating in the galactic halos move with a 
velocity of the order vjc ~ 10""^, hence these outgoing photons (as any particle in any of the 
allowed two body final states) will then be nearly monochromatic, with energy of the order 
of the neutralino mass |llll 11121 11181 11021 IIOOL l(i()| . There is no other known astrophysi- 
cal source with such a signature: the detection of a line signal out of a spectrally smooth 
gamma-ray background would be a spectacular confirmation of the existence of dark matter 
in form of exotic massive particles. 

6.5.1 Sources and fluxes 

Following the discussion in |(i()| . the monochromatic gamma-ray flux measured in a detector 
with angular acceptance Af2 is: 

$,(^, A^) . 0.94 . 10-" ( ^oXm^-O (^) ^ ^ ' ^^'^^^'' ' ^"^ '^"" '" "" ' 

(41) 
where ip is an angle to label the direction of observation and where iV^ — 2 for X X ^ 7 7, 
iV^ = 1 for XX ^ ■2^7- Here the dimensionless line-of-sight dependent function is 

1 /- 1 ^^ 



and its angular average over the resolution solid angle Afi is 

{JmAn^^J^jn'Ji^,') , (43) 

Analogously, the gamma-ray flux with continuum energy spectrum is obtained by replac- 
ing Nj vaxo-y with ^ , dNji /dE vaj, where the sum is over all tree level final states. 

Finally, the formalism we introduced can be used also to estimate the fiux in the simple 
case of a single source which, for a given detector, can be approximated as point-like. If such 
a source is in the direction -0 at a distance d, Eq. H43|l becomes: 

<^(^)>- - 8:51;^ ■ ( o.3GeV/cm3 )^ ' i ' ^ j ''^ '^^^ ^''^ 

where the integral is over the extension of the source (much smaller than d). 
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Several targets have been discussed as sources of gamma-rays from the annihilation of 
dark matter particles. An obvious source is the dark halo of our own galaxy |114| and in 
particular the Galactic center, as the dark matter density profile is expected, in most models, 
to be peaked towards it, possibly with huge enhancements close to the central black hole. The 
Galactic center is an ideal target for both ground- and space-based gamma-ray telescopes. 
As satellite experiments have much wider field of views and will provide a full sky coverage, 
they will test the hypothesis of gamma-rays emitted in clumps of dark matter which may 
be present in the halo |115l II 101 lllfil I117J . Another possibility which has been considered is 
the case of gamma-ray fluxes from external nearby galaxies |118| . Furthermore, it has been 
proposed to search for an extragalactic flux originated by all cosmological annihilations of 
dark matter particles [1191 I12UJ . 

DarkSUSY is suitable to compute the gamma-ray flux from all these (and possibly other) 
sources. Two cases are fully included in the package: assuming that neutralinos are smoothly 
distributed in the Galactic halo with p^ equal to the dark matter density profile, in DarkSUSY 
Eq. (|43|l is computed for a specified halo profile and any given ij) and Af2 [HO]. The second 
option deals with the case of a portion of dark matter being in the form of clumps, each of 
which is treated as a non-resolvable source in the detector, distributed in the Galaxy according 
to a probability distribution which can be specified by the user; in DarkSUSY the default 
probability distribution stems from the assumption that it follows the dark matter density 
profile (see |llfi| for details). It is straightforward to extend this to all other astrophysical 
sources; in case of cosmological sources one has just to pay attention to include rcdshift 
effects and absorption on starlight and infrared background, see |120| . 

The case of the possible enhancement, a "spike" in the vicinity of the galactic center |121| 
should be kept in mind. However, since there is no consensus in the literature |122l 11281 I124J 
about important quantities for the calculation such as the magnetic field radial profile and 
the optical depth for synchrotron self-absorption, we have chosen not to include routines for 
the effects of this possible enhancement of gamma rays and neutrinos. And the very existence 
of a spike is dependent on fine details, still unknown for the Milky Way [12511126] . 

6.5.2 Gamma rays with continuum energy spectrum 

The gamma-ray flux produced in neutralino annihilations through tt*^ decays can be large but 
in general lacks distinctive features. The photon spectrum in the process of a pion decaying 
into 27 is, independent of the pion energy, peaked at half of the 7r° mass, about 70 MeV, 
and symmetric with respect to this peak if plotted in logarithmic variables (i.e. dNj/dE vs. 
log_E). Of course, this is true both for pious produced in neutralino annihilations and, e.g., 
for those generated by cosmic ray protons interacting with the interstellar medium. 

When considered together with the cosmic ray induced Galactic gamma-ray background, 
the neutralino induced signal looks like a component analogous to the secondary flux due to 
nucleon nucleon interactions: it is dwarfed by the bremsstrahlung component at low energy, 
while it may be the dominant contribution at energies above 1 GeV or so. In DarkSUSY 
the continuum gamma flux from all annihilation channels is computed and may be easily 
obtained for a given energy or energy threshold. 

6.5.3 XX -^ 77 

At the one-loop level, it is possible to get two-body final states containing one or two photons, 
with a distinctive signature which may provide a "smoking gun" for dark matter. The 
drawback is of course that the processes are loop-suppressed, so one probably needs a halo 
with a large central concentration, or small-scale structure ("clumps" of dark matter) to 
detect a signal. 
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In DarkSUSY the full expression for the annihilation cross section of the process 

X? + X? ^ 7 + 7 (45) 

is computed at full one loop level, in the limit of vanishing relative velocity of the neutralino 
pair, i.e. the case of interest for neutralinos in galactic halos; the outgoing photons have an 
energy equal to the mass of Xi- 

E^^m^. (46) 

The total amplitude is implemented in DarkSUSY as the sum of the contributions obtained 
from four different classes of diagrams: 

A = Aff + Ah+ +Aw+ Ag, 

where the indices label the particles in the internal loops, i.e., respectively, fermions and 
sfermions, charged Higgs and charginos, W-bosons and charginos, and, in the gauge we chose, 
charginos and Goldstone bosons. For every A term, real and imaginary parts are separately 
computed; the full set of analytic formulas are given in 102^, following the notation of jli;-{| . 
where some of the contributions were first computed. They are rather lengthy expressions 
with non trivial dependences on various combinations of parameters in the MSSM. 

The branching ratio for neutralino annihilations into 27 is typically not larger than 1%, 
with the largest values of va2-y , for neutralinos with a cosmologically significant relic abun- 
dance, in the range 10~^^-10~^^ cm^s~^. Such values may be large enough for the discovery 
of this signal in upcoming measurements; at the same time it should be kept in mind that 
very low values for the cross section are possible as well. 

In the very high mass range above a TeV, it has been suggested that the line rates may 
be very much larger due to binding effects and resonant conversion between neutralinos and 
charginos |127j . In the present version of DarkSUSY we have not included these effects, 
however. 

6.5.4 XX -" Z-f 

The process of neutralino annihilation into a photon and a Z° boson |1U(JJ 

X? + X?^7 + ^° (47) 

also gives a nearly monochromatic line (with a small smearing caused by the finite width 

of the Z'^ boson, in any case unlikely to be important for current or proposed gamma ray 

experiments), with energy 

2 

£; =M --^. (48) 

The steps followed in DarkSUSY to compute the cross section are essentially the same 
as those described for the 27 case. Again the total amplitude is obtained by summing the 
contribution from four classes of diagrams and by splitting for each of them real and imaginary 
parts. The analytic formulas were derived in '100', and are much less compact than those 
obtained for the process of neutralino annihilation into two photons. 

The maximum value of vaz-y , for neutralinos with a cosmologically significant relic abun- 
dance, is around 2 • 10~^^ cm'^s"^ and occurs for nearly pure higgsinos. In the heavy mass 
range, the value of vaz-f reaches a plateau of around 0.6 • 10~^^ cm^s^^. This interesting 
effect of a non-diminishing cross section with Higgsino mass (which is due to a contribution 
to the real part of the amplitude) is also valid for the 27 final state in the corresponding 
limit, with a value of 1 • 10^^* cm^s~^ |1U2| . Since the gamma-ray background drops rapidly 
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with increasing photon energy, these processes may be interesting for detecting dark matter 
neutralinos near the upper range of aUowed neutrahno masses. 

Whenever the hghtest neutrahno contains a significant Wino or Higgsino component the 
value of vaz-y may be as large as, or even larger than, twice the value of va2'y It is therefore 
usually not a good approximation to neglect the Zj state compared to 27. 

6.6 Neutrinos from the halo 

Usually, the flux of neutrinos from annihilation of neutralinos in the Milky Way halo is 
too small to be detectable, but for some clumpy or cuspy models, or for annihilation in a 
possible spike around the central black hole, it might be detectable. The calculation of the 
neutrino-flux follows closely the calculation of the continuous gamma ray flux, with the main 
addition that neutrino interactions close to the detector are also included. Hence, both the 
neutrino flux and the neutrino-induced muon flux can be obtained. The neutrino to muon 
conversion rate in the Earth can also be obtained. The neutrino rate from other sources than 
the interior of the Earth or the Sun is generally negligible. If there would exist a spike at the 
galactic center |121| . there may in some cases be a detectable flux. These neutrino rates are 
constrained by existing limits on the gamma-ray flux J121I I128J . 

7 Examples of results obtained with DarkSUSY 

We will here go through a set of benchmark models as examples of the performance of 
DarkSUSY. We will consider two popular setups. We will start with a set of mSUGRA models 
and then turn to more general MSSM models. We will here use the default DarkSUSY setup, 
in particularan isothermal sphere with a Maxwell-Boltzmann velocity distribution for the 
halo model. 

7.1 Benchmark models in mSUGRA 

We will consider here some of the benchmark models from Battaglia et al. |129j . In table ^ 
we list the properties of the selected models, as derived by DarkSUSY using to the ISASUGRA 
7 . 69 package to describe the renormalization group running of the theory from the grand 
unification scale to the low energy scale. The models we are focusing on are those with a top 
mass of 175 GeV and that are stiU viable with ISASUGRA 7.69. As the table in |129| was 
produced with ISASUGRA 7.67, some differences occur due to different versions of the codes, 
e.g. model M is no longer physical in ISASUGRA 7.69, and thus not included here. 

Our results for the relic density and for the SUSY contributions to B{b — + 57) and to the 
anomalous magnetic moment of the muon a^, agree reasonably well with those in 129 . The 
expected sensitivities of future neutrino telescopes are of the order of 20 events km~^ yr~^ for 
the Earth and 50-1000 events km~^ yr""'^ for the Sun; we see in the table that all benchmark 
models, unfortunately, produce much lower fluxes. The outreach of future direct detection 
experiments is expected to be of the order of 10~^ pb for the spin-independent scattering 
cross section; we see that some of models we are considering are potentially detectable. The 
cross section for annihilation into gamma rays (times the number of photons produced) are 
also given for these models; the detectability depends strongly on the halo profile, but one can 
in general say that these cross sections are too low to be seen in current data, unless the halo 
profile is very cuspy towards the galactic center. The e"*" fluxes are here given in the energy bin 
6.0-8.9 GeV of the HEAT 94-1-95 130; experiment. The measured e+ flux is (7.2±1.2) x 10^6 
GeV~^ cm~^ s~^ sr~^, i.e. the predicted fluxes are much lower than the measured one. For 
the antiprotons, we choose, as an example, to show the predicted (average) flux in the energy 
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Table 1: Relic density and various rates for the benchmark models of 129. . There are five free parameters in mSUGRA: tan/?, sign(/x), 
mi/2, rno and Aq. The latter three are the unification values (at the grand unification scale) of the soft supersymmetry breaking fermionic 
mass parameters, scalar mass parameters and trilinear scalar coupling parameters, respectively. All the benchmark models have ^0 = 0- 
We have here used ISASUGRA 7.69 for the RGE-running (but have not taken the b, t, and r Yukawa couplings from ISASUGRA). The 
neutrino-induced muon fluxes from the Earth and the Sun are for a threshold of 1 GeV. 'best' refers to the suppressed fluxes resulting from 
the estimate in |H1] and 'gauss' refers to the usual approximation of a Gaussian velocity distribution of neutralinos for capture in the Earth. 
The scattering cross sections are given here evaluated both with the standard expressions (labelled 'std') and with the Drees and Nojiri 
expressions |70| (labelled 'dn'). The continuum 7's are given in terms of the number of photons times the cross section, N cont ('''^)' ^^'^ 
refers to 7's above 1 GeV. The e+ flux is the average solar modulated flux in the energy range 6.0-8.9 GeV, i.e. in one of the HEAT 94-95's 
binsll30l. The p flux is the average solar modulated flux in the energy range 0.56-0.78 GeV, i.e. in one of the BESS 98 bins |131| . The 
D flux is the average flux in the energy range 0.1-0.4 GeV, as applicable to e.g. the proposed GAPS probe \IIV2\ . In this case, the flux is 
the average of the solar modulated fluxes at solar minimum and maximum. Parameters other than those mentioned above are set to their 
default values as described in the text. 



range 0.56-0.78 GeV, which is one of the BESS 98 |131j bins. The measured flux in this 
bin is (1.231q'33) ^ 10~® GeV^^ cm^^ s^^ sr^^, i.e. the measured flux is much higher than 
the predicted one in this energy range. For the antideuterons, we show the expected flux 
in the energy range 0.1-0.4 GeV, which is reasonable for the proposed GAPS probe 132]. 
For the antideuterons, there is essentially no background and the sensitivity is thus given by 
the ability to detect one antideuteron. For GAPS this corresponds to a sensitivity of about 
2.6 X 10"^'^ GeV~^ cm~^ s~^ sr~^. Hence, two of these benchmark models have high enough 
fluxes for being just about detectable in this way. For further examples of rates in mSUGRA, 
as calculated with DarkSUSY , see e.g. |133| . 

7.2 Benchmark models in MSSM 

We will now turn to models in the MSSM framework as generated by fixing free parameters 
at the weak scale. We refer to the setup with seven free parameters, i.e. /i, M2, niA, tan/3. 
Too, At and Ab, we described in Section 12 We will here show an example of results that 
can be obtained with a simple scan over this parameter space. We have generated 5000 
models assuming: ^l G [-3000,3000] GeV, Ah G [-3000,3000] GeV, toa G [100,1000] GeV, 
tan/3 G [1,55], toq G [100,5000] GeV, At G [-3,3]too and Ab G [-3,3]too. We have then 
selected a few sample models with a relic density in the range 0.09 < Jl^/i^ < 0.11 and with 
reasonably high detection rates; these are the first ten models in Table|21 In addition to this 
scan, we have also made a small scan of 300 models in which we required the mass of the 
CP-odd Higgs boson. A, to be in the range toa G [90, 150] GeV; model number 11 in Table 
121 has been chosen from this latter scan to illustrate that higher-rate models are possible to 
find with these kind of dedicated scans. 

As seen in the table, the rates are typically slightly higher than for the mSUGRA bench- 
mark models considered above. Please remember though, that this set of MSSM benchmark 
models is achieved with a rather small scan over the MSSM parameter space. Models with 
even higher rates would be possible to find with more extensive scans of the parameter space. 

Many of our models produce fiuxes from the Sun that exceed the future sensitivity of 
50-1000 events km~^ yr~^ and hence would be detectable. For the Earth, on the other hand, 
the fluxes in this set of models are typically much lower than the projected future limit of 20 
events km~^ yr~^. Most of the selected models are potentially detectable with future direct 
detection experiments as they have a spin-independent scattering cross section larger than 
10^^ pb. Note the complementarity between the direct detection signal and the neutrino 
flux from the Sun. For example, for model 9, the spin-independent scattering cross section 
is very close to the sensitivity of future detectors, whereas the neutrino-flux from the Sun 
is clearly above expected future sensitivities of neutrino telescopes like Antares or IceCube. 
The gamma-ray yields are in general larger than for the mSUGRA benchmark points, but the 
issue regarding detectability is still difficult to address as it is strongly affected by the halo 
model dependence. The positron and antiproton fluxes are generally enhanced as well, but 
still well below measured fluxes. Finally, there are several models among those we selected 
which have an antideuteron flux in the energy range 0.1-0.4 GeV exceeding 2.6 x 10~^^ 
GeV^^ cm^^ s^^ sr~^, the expected sensitivity of the GAPS detector |132| . 

7.3 Discussion 

The potential of DarkSUSY package has been illustrated here, for a few sample models and 
in two popular scenarios. The code provides the state of the art calculation of the neutralino 
relic abundance, and allows for the estimate of several quantities of interest for dark matter 
searches. It has been shown that the code has a very flexible structure in the definition of 
the supersymmetric dark matter candidate, as well as a rather broad freedom in the choice 
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8.33 ■ lO-'^ 

3.68 ■ 10-^ 

1.57 ■ 10 



4.20 
-0.24 

■^2 ^1,2,3 

0.1032 
1.62 ■ 10-*^ 
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5.26 

3.08 

3.22 

5.26 

3.22 

5.26 

3.45 ■ 1 

5.07 

4.91 

1.91 

8.09 ■ 10 

2.06 ■ 10 
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■ 10"" 
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■ 10"" 
10" 



10" 



10" 



4.51 

0.50 

Y+ y" 
^2 ^1,2,3 

0.0977 

4.67 ■ 10"^ 

-4 



9.79 ■ 10 

1.74 ■ 10 
1.01 ■ 10" 
6.10 ■ 10" 

1.01 ■ 10" 

6.02 ■ 10" 
8.04 ■ 10' 

0.0077 

5.41 

4.20 ■ 10"" 

1.25 ■ 10"" 

6.55 ■ 10- 



2.82 

15.87 

Y+ y" 

0.0968 

2.17 ■ 10-'- 



2.04 

3.66 

4.14 ■ 

3.47 ■ 

4.14 ■ 

3.40 ■ 

4.64 

0.059 

0.11 

7.78 ■ 10 

2.05 ■ 10 

4.33 ■ 10" 
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Table 2: Same as Tabled but for the MSSM benchmark models of our sample scan. 



of the relevant astrophysical setups. Outputs are in simple and general formats, which, as 
we have shown, can be very easily compared to current and future sensitivities. Some trends 
on predictions can be extracted from benchmark models, as we have done for some of the 
mSUGRA models proposed in |129j , and with sample models selected according to their relic 
abundance in a more general low energy scan (with the latter being more promising than the 
former) . One should keep in mind however that firmer statements are possible just in light 
of dedicated and more extensive scans. 

8 Conclusions 

We have here described the computer package DarkSUSY, that can be used to calculate various 
quantities of interest for supersymmetric dark matter searches. We have gone through great 
efforts and used state-of-the-art techniques to obtain a package that can deliver very accurate 
results in a flexible setting. We also believe that this package can be of great use for the 
physics community. 

In this paper we have described briefly the ingredients of DarkSUSY and shown, with some 
examples, what one can calculate with it. We encourage the reader to download DarkSUSY 
and start using it. 
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A Included coannihilations 

In this Appendix, we tabulate the coannihilation processes computed in DarkSUSY. For more 
details, see jMi- 

In the table, we list all 2 — > 2 tree-level coannihilation processes with sfermions, charginos 
and neutralinos. All the processes are included in the DarkSUSY code. 

It should be noted that we have not included all flavour-changing charged current dia- 
grams. The DarkSUSY vertex code for the charged current couplings is written in a gen- 
eral form that includes all possible flavour-changing (and flavour-conserving) vertices. The 
flavour-conserving couplings are much larger than the flavour-changing. For the sfermion 
coannihilations with charged currents we only take the flavour-conserving contributions, while 
for the chargino coannihilations we include the flavour-changing contributions as well. In a 
future version of DarkSUSY, we may as well include the flavour-changing processes for the 
sfermion coannihilations, even if they are not expected to be important. 

We have used the notation / for sfermions and / for fermions. Whenever the isospin of 
the sfermion/fermion is important, it is indicated by an index u (Ta = 1/2) or d (T^ — —1/2). 
The sfermions have an additional mass eigenstate index, that can take the values 1 and 2 
(except for the sneutrinos which only have one mass eigenstate). A further complication 
to the notation is when the sfermions and fermions in initial, flnal and exchange state can 
belong to different families. Primes will be used to indicate when we have this freedom to 
choose the flavour. So, e.g. /„ and /„ will belong to the same family while /„ and /^ can 
belong to the same or to different families. Note that the colour index of (s)quarks as well 
as gluons (g) and gluinos {g) is suppressed. 
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running from 1 to 4. The charginos are similarly denoted Xj with the index taking the values 



Besides the sfermions we also have neutralinos and charginos in the initial states. The 
notation used for these are the following. The neutralinos are denoted by x'j with the index 

from 1 to 4. The charginos are similarly denoted y^ 
1 and 2. 

In the table, a common notation is introduced for gauge and Higgs bosons in the final 
state. We denote these with B with an upper index indicating the electric charge. So B^ 
means i?°, H2, H^, Z, 7 and g while B^ is H^ and W^ . We will use additional lower indices 
m and n when we have more than one boson in the final state. Thus indicating that the 
bosons can be either different or identical. Note that the case of two different bosons also 
includes final states with one gauge boson and one Higgs boson. 

The table has been made very general. This means that when a set of initial and final 
state (s)particles has been specified, the given process might not run through all the exchange 
channels listed for the generic process. Exceptions occur whenever an exchange (s)particle 
does not couple to the specific choice of initial and/or final state. As an example we see that 
since the photon does not couple to neutral (s)particles, none of the exchange channels listed 
for the generic process /i + X? -^ B^ + f actually exist for the specific process v + x^ ~^ j + i^- 
All these exceptions can be found in the extended tables in Ref. |184| . Also note that the 
list of processes is not complete with respect to trivial charge conjugation. For each process 
of non-vanishing total electric charge in the initial state there exists another process which 
is obtained by charge conjugation. 
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Diagrams 
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Table 3: Included coannihilation processes through s— , t— , m— channels and four-point in- 
teractions (p). For the fdifd* processes the corres 
be obtained by interchanging the u and d indices. 



teractions (p). For the fdifd* processes the corresponding process for up- type sfermions can 
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